# TODO: create PM image from aod and temp (PM regression)
# 
# Author: phamha
###############################################################################

#import library
library(gstat)
library(base)
library(RPostgreSQL)
library(stringr)
library(raster)
library(gdalUtils)
library(rgdal)
library(rPython)

host_name = '192.168.0.4'
database_name = 'fimo'
user_name = 'rasdaman'
password = 'rasdaman'

min_aod_mod = 0.0015
max_aod_mod = 1.165031427
min_temp_mod = 291.8899935
max_temp_mod = 317.5399929

min_aod_myd = 0.057872928
max_aod_myd = 1.509154063
min_temp_myd = 281.2699937
max_temp_myd = 316.6299929


min_aod_npp = 0.065608211
max_aod_npp = 1.801370621
min_temp_npp = 281.0396729
max_temp_npp = 326.4881287


data_folder = "/var/www/html/"
met_folder = "/home/phamha/MetNPP/"
myd_folder = "/var/www/html/fimo/apom/Product/MYD01PM/"
mod_folder = "/var/www/html/fimo/apom/Product/MOD01PM/"
shape_file = "/home/phamha/BDNEN/VNM_adm0.shp"
#shape_file = "E:/FIMO/DATA/BDNEN/VNM_adm0.shp"

mod_time_query = "select aqstime from res.satresampviirs where filename like '%"

mod_query = "SELECT mod04.aqstime, mod04.filename,mod04.filepath,mod07.filename as temp_filename,mod07.filepath as temp_filepath FROM res.satresampviirs as mod04 inner join res.satresampviirstemperature as mod07 ON (mod04.aqstime = mod07.aqstime) where mod04.aqstime < '2014-12-31 00:00:00'::timestamp"


regress_predict = function(sate_data,model_type,cook_type,aod,temp,avg_temp,avg_rh,avg_preci){
	if(sate_data=="mod"&&model_type=="linear"&&cook_type=="4np"){
		pm25 = 22.58132724* aod + 0.08452807 * temp + (-2.98371028)*avg_temp + (-0.58041657)*avg_rh + 0.64525692 * avg_preci + 113.98273058
	}
	if(sate_data=="mod"&&model_type=="linear"&&cook_type=="fdist"){
		pm25 = 29.3009026972483* aod - 0.0121947942915479 * temp + (-4.05178592218209)*avg_temp + (-1.17578020320201)*avg_rh + 1.01861037054962 * avg_preci + 213.045850613292
	}
	if(sate_data=="myd"&&model_type=="linear"&&cook_type=="4np"){
		pm25 = 25.6965586 * aod + 0.1525928 * temp + (-2.1894375)*avg_temp + (-0.6837351)*avg_rh + 0.4189560 * avg_preci + 70.1740458
	}
	if(sate_data=="myd"&&model_type=="linear"&&cook_type=="fdist"){
		pm25 = 30.5744605532908 * aod + 0.118323308557864 * temp + (-2.47362345170908)*avg_temp + (-0.795065245241148)*avg_rh + 0.45323741215544 * avg_preci + 94.9758204191974
	}
	if(sate_data=="npp"&&model_type=="linear"&&cook_type=="4np"){
		pm25 = 14.9045673 * aod + 0.0647682902577907 * temp + (-2.91326445481888)*avg_temp + (-1.20281335981164)*avg_rh + 0.446332747336052 * avg_preci + 154.899079958179
	}
	
	return(pm25)
}
#get data from DB
getDataFromDB = function(sql_command){
	driver = dbDriver("PostgreSQL")
	connect = dbConnect(driver,dbname = database_name,host = host_name,port=5432,user = user_name,password= password)
	rs = dbSendQuery(connect,sql_command)
	data=fetch(rs,n=-1)
	return (data)
}

create_regress_pm = function(sate_data,model_type,cook_type,aod_file,temp_file){
	if(sate_data=="mod"){
		time_query = mod_time_query
		min_aod = min_aod_mod
		max_aod = max_aod_mod
		min_temp = min_temp_mod
		max_temp = max_temp_mod
		
	}
	if(sat_data=="myd"){
		time_query = myd_time_query
		min_aod = min_aod_myd
		max_aod = max_aod_myd
		min_temp = min_temp_myd
		max_temp = max_temp_myd
	}
	if(sat_data=="npp"){
		time_query = npp_time_query
		min_aod = min_aod_npp
		max_aod = max_aod_npp
		min_temp = min_temp_npp
		max_temp = max_temp_npp
	}
	
	start_index = regexpr("M[^M]*$",aod_file)
	end_index = regexpr(".tif",aod_file) - 1
	file_name = substr(aod_file,start_index,end_index)
	
	time_query = paste(time_query,file_name,"%'",sep="")
	data = getDataFromDB(time_query)
	
	mod04_aqstime = data$aqstime[1]
	aqstime = strptime(mod04_aqstime,format="%Y-%m-%d %H:%M:%S")
	aqstime = aqstime + 25200	
	month = format.Date(aqstime,"%m")
	year = format.Date(aqstime,"%Y")

	
	shape_dataset <- system.file(shape_file, package="gdalUtils")
	aod_file = "C:/MOD04_L2.A2014024.0310.051.2014028020029.hdf_DT_10km.tif"
	gdal_setInstallation(ignore.full_scan=FALSE)
	gdal_rasterize(shape_dataset,aod_file,b=1,i=TRUE,burn=-9999,l="VNM_adm0")
	
	aod_dataset = raster(aod_file)
	aod_data = values(aod_dataset)
	#aod_data = aod_data * 0.00100000004749745
		
	corxy = coordinates(aod_dataset)
	x = corxy[,'x']
	y = corxy[,'y']
		

	temp_dataset = raster(temp_file)
	temp_data = values(temp_dataset)
	#temp_data = (temp_data + 15000) * 0.00999999977648258
		
	avg_temp_file  = paste(met_folder,"temp",as.numeric(month),".tif",sep="")
	avg_rh_file    = paste(met_folder,"rh",as.numeric(month),".tif",sep="")
	avg_preci_file = paste(met_folder,"preci",as.numeric(month),".tif",sep="")
		
		
	avg_temp_dataset  =  raster(avg_temp_file)
	avg_rh_dataset    =  raster(avg_rh_file)
	avg_preci_dataset =  raster(avg_preci_file)
		
	avg_temp_data  = values(avg_temp_dataset) 
	avg_rh_data    = values(avg_rh_dataset) 
	avg_preci_data = values(avg_preci_dataset)
		

	pm25 = regress_predict(sate_data,model_type,cook_type,aod_data,temp_data,avg_temp_data,avg_rh_data,avg_preci_data)

	table = data.frame(x,y,aod_data,temp_data,avg_temp_data,avg_rh_data,avg_preci_data,pm25)
	table$pm25[table$aod_data<min_aod|table$aod_data>max_aod|table$temp_data<min_temp|table$temp_data>max_temp]<-NA
	og_raster = aod_dataset
	totalCell = ncell(og_raster)
	og_raster[1:totalCell] = table$pm25

	pm_file = str_replace(aod_file,".tif","_.tif")
	writeRaster(og_raster,filename=pm_file,format="GTiff")
	print(file_name)
			
}
createKrigingImage = function(regressPm_file){
	
	#PM values
	pmRaster=raster(regressPm_file)
	pm=values(pmRaster)
	corxy=coordinates(pmRaster)
	x=corxy[,'x']
	y=corxy[,'y']
	
	
	totalCell=length(pmRaster)
	cell = c(1:totalCell)
	
	table=data.frame(cell,x,y,pm)
	newTable=table
	testTable=subset(table,pm<0)
	trainTable=subset(table,pm>=0)
	
	#caculate variogram
	empiVario=variogram(pm~1,locations=~x+y,data=trainTable)
	
	#sph fit
	sill=min(empiVario$gamma)
	#sphModel=vgm(psill=sill,model="Sph",nugget=0,range=min(empiVario$dist))
	sphModel=vgm(model="Sph",nugget=0,range=1)		
	sphFit=fit.variogram(empiVario,sphModel)
	
	universal_result=krige(id="pm",formula=pm~x+y,data=trainTable,newdata=newTable,model=sphFit,locations=~x+y)
	
	#edit tiff
	universalPMRaster=pmRaster
	universalPMValue=universal_result[,3]
	universalPMRaster[1:totalCell]=universalPMValue
	
	#edit error tiff
	errorPMRaster=pmRaster
	errorPMValue=universal_result[,4]
	errorPMRaster[1:totalCell]=errorPMValue
	
	#save uk result to tiff
	uk_file = str_replace(regressPm_file,"rg.tif","uk.tif")
	writeRaster(universalPMRaster,filename=uk_file,format="GTiff")
	
	#save uk error to tiff
	error_file = str_replace(regressPm_file,"rg.tif","error.tif")
	writeRaster(errorPMRaster,filename=error_file,format="GTiff")
	print(uk_file)
	
}
aod_file = "C:/MOD04_L2.A2014024.0310.051.2014028020029_rg.tif"
temp_file = "C:/MOD04_L2.A2014024.0310.051.2014028020029_rg.tif"
create_regress_pm("mod","linear","4np",aod_file,temp_file)
print("done")
